home *** CD-ROM | disk | FTP | other *** search
/ Packard Bell - Internet on a CD / internet on a cd.cdr / Internet / sites / Clementine_NASA / clemdsrc.hqx / decomp.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-06-27  |  19.3 KB  |  738 lines

  1. /*
  2.  *    THIS ROUTINE IS PART OF THE CLEMENTINE PDS FILE READER PROGRAM.
  3.  *    IT WAS WRITTEN BY ACT CORP. IN DIRECT SUPPORT TO THE 
  4.  *    CLEMENTINE (DSPSE) PROGRAM.
  5.  *
  6.  *    IF YOU FIND A PROBLEM OR MAKE ANY CHANGES TO THIS CODE PLEASE CONTACT
  7.  *    Dr. Erick Malaret at ACT Corp. 
  8.  *            tel: (703) 742-0294
  9.  *                 (703) 683-7431
  10.  *                       email:    nrlvax.nrl.navy.mil
  11.  *    
  12.  *
  13.  */
  14. /*
  15.     Program to calculate the Discrete Cosine Transform of a 8x8 block of data
  16. */
  17.  
  18. #include <stdio.h>
  19. #ifdef __TURBOC__ 
  20. #include <alloc.h>
  21. #else
  22. #include <stdlib.h>
  23. #endif
  24. #include <string.h>
  25. #include <math.h>
  26. #include "jpeg_c.h"
  27.  
  28. #define TRUE 1
  29. #define FALSE 0
  30.  
  31. float    qtable[64];
  32. float   U[64];
  33. short   ULS[64];
  34. int    zzseq[] = {0,1,8,16,9,2,3,10,17,24,32,25,18,11,4,5,12,19,26,33,40,
  35.            48,41,34,27,20,13,6,7,14,21,28,35,42,49,56,57,50,43,36,29,
  36.            22,15,23,30,37,44,51,58,59,52,45,38,31,39,46,53,60,61,54,
  37.            47,55,62,63};
  38.  
  39. void core(void);
  40.  
  41.  
  42. void decomp(BitStream *bs,CHARH *Image,long rows,long cols)
  43. {
  44.     short   i,j,k,indexi,indexip8,indexj,indexjp8;
  45.     float   temp;
  46.     long    rowsleft, colsleft, icols;
  47.     short   Done, Diff, Pred;
  48.  
  49.     /*********************** Image Decompression *********************/
  50.     Done = FALSE;
  51.     rowsleft = rows;
  52.     colsleft = cols;
  53.     indexj = 0;
  54.     indexi = 0;
  55.     Pred = 0;
  56.     
  57.     while ( !Done ) {
  58.         /* Decode block data */
  59.         decode(ULS,bs);
  60.         Diff = ULS[0] + Pred;
  61.         Pred = Diff;
  62.         ULS[0] = Diff;
  63.  
  64.         /* Calculate Inverse Discrete Cosine Transform */
  65.         /*dequantize(U,ULS1,pt_qtable);  64 multiplications */
  66.         for (i=0; i<64; i++) U[zzseq[i]] = ULS[i] * qtable[i];
  67.         
  68.         core();
  69.  
  70.         /*ilevelshift(U);*/
  71.         for (i=0; i<64; i++) {
  72.                         U[i] += 128.0;
  73.             U[i] = floor( U[i] + 0.5 );
  74.             if ( U[i] > 255.0 ) U[i] = 255.0;
  75.             else if ( U[i] < 0.0 ) U[i] = 0.0;
  76.             else;
  77.         }
  78.         
  79.         indexip8 = indexi + 8;
  80.         indexjp8 = indexj + 8;
  81.  
  82.         if ( (rowsleft > 8) && (colsleft > 8) ) {
  83.  
  84.             for (i=indexi, k=0; i < indexip8; i++)
  85.                 for (j=indexj, icols=i*cols; j < indexjp8; j++, k++)
  86.                     Image[icols+j] = (char)U[k];
  87.  
  88.             indexj += 8;
  89.             colsleft -= 8;
  90.         }
  91.         else if ( (rowsleft > 8) && (colsleft <= 8) ) {
  92.  
  93.             for (i=indexi, k=0; i < indexip8; i++) {
  94.                 for (j=indexj, icols=i*cols; j < cols; j++, k++)
  95.                     Image[icols+j] = (char)U[k];
  96.                 k += (8 - colsleft);
  97.             }
  98.  
  99.             indexj = 0;
  100.             indexi += 8;
  101.             rowsleft -= 8;
  102.             colsleft = cols;
  103.         }
  104.         else if ( (rowsleft <= 8) && (colsleft > 8) ) {
  105.  
  106.             for (i=indexi, k=0; i < rows; i++)
  107.                 for (j=indexj, icols=i*cols; j < indexjp8; j++, k++)
  108.                     Image[icols+j] = (char)U[k];
  109.  
  110.             indexj += 8;
  111.             colsleft -= 8;
  112.         }
  113.         else {
  114.  
  115.             for (i=indexi, k=0; i < rows; i++) {
  116.                 for (j=indexj, icols=i*cols; j < cols; j++, k++)
  117.                     Image[icols+j] = (char)U[k];
  118.                 k += (8 - colsleft);
  119.             }
  120.  
  121.             Done = TRUE;
  122.         }
  123.     }
  124. }
  125.  
  126.  
  127. void core(void)
  128. {
  129.     short   i;
  130.     float   out[64], out1[64];
  131.     float   temp1,temp2,temp3,temp4,temp5,temp6,temp7,temp8;
  132.     float   dummy1, dummy2, dummy3;
  133.     float   buffer11,buffer12,buffer13,buffer14,buffer15,buffer16,buffer17,buffer18;
  134.     float   buffer21,buffer22,buffer23,buffer24,buffer25,buffer26,buffer27,buffer28;
  135.  
  136.     /******** Preadditions ********/
  137.     out[0] = U[0];
  138.     out[1] = U[32];
  139.     out[2] = U[16] - U[48];
  140.     out[3] = U[16] + U[48];
  141.     dummy1 = U[8] - U[56];
  142.     dummy2 = U[24] - U[40];
  143.     out[4] = dummy1 - dummy2;
  144.     out[5] = dummy1 + dummy2;
  145.     out[6] = -U[8] - U[56];
  146.     out[7] = U[24] + U[40]; /* 8 additions */
  147.     
  148.     out[8] = U[4];
  149.     out[9] = U[36];
  150.     out[10] = U[20] - U[52];
  151.     out[11] = U[20] + U[52];
  152.     dummy1 = U[12] - U[60];
  153.     dummy2 = U[28] - U[44];
  154.     out[12] = dummy1 - dummy2;
  155.     out[13] = dummy1 + dummy2;
  156.     out[14] = -U[12] - U[60];
  157.     out[15] = U[28] + U[44]; /* 8 additions */
  158.  
  159.     temp1 = U[2] - U[6];
  160.     temp2 = U[34] - U[38];
  161.     temp3 = U[18] - U[22];
  162.     temp4 = U[50] - U[54];
  163.     temp5 = U[10] - U[14];
  164.     temp6 = U[26] - U[30];
  165.     temp7 = U[58] - U[62];
  166.     temp8 = U[42] - U[46];
  167.     out[16] = temp1;
  168.     out[17] = temp2;
  169.     out[18] = temp3 - temp4;
  170.     out[19] = temp3 + temp4;
  171.     dummy1 = temp5 - temp7;
  172.     dummy2 = temp6 - temp8;
  173.     out[20] = dummy1 - dummy2;
  174.     out[21] = dummy1 + dummy2;
  175.     out[22] = -temp5 - temp7;
  176.     out[23] = temp6 + temp8; /* 16 additions */
  177.  
  178.     temp1 = U[2] + U[6];
  179.     temp2 = U[34] + U[38];
  180.     temp3 = U[18] + U[22];
  181.     temp4 = U[50] + U[54];
  182.     temp5 = U[10] + U[14];
  183.     temp6 = U[26] + U[30];
  184.     temp7 = U[58] + U[62];
  185.     temp8 = U[42] + U[46];
  186.     out[24] = temp1;
  187.     out[25] = temp2;
  188.     out[26] = temp3 - temp4;
  189.     out[27] = temp3 + temp4;
  190.     dummy1 = temp5 - temp7;
  191.     dummy2 = temp6 - temp8;
  192.     out[28] = dummy1 - dummy2;
  193.     out[29] = dummy1 + dummy2;
  194.     out[30] = -temp5 - temp7;
  195.     out[31] = temp6 + temp8; /* 16 additions */
  196.  
  197.     buffer11 = U[1] - U[7]; buffer21 = U[3] - U[5];
  198.     buffer12 = U[33] - U[39]; buffer22 = U[35] - U[37];
  199.     buffer13 = U[17] - U[23]; buffer23 = U[19] - U[21];
  200.     buffer14 = U[49] - U[55]; buffer24 = U[51] - U[53];
  201.     buffer15 = U[9] - U[15]; buffer25 = U[11] - U[13];
  202.     buffer16 = U[25] - U[31]; buffer26 = U[27] - U[29];
  203.     buffer17 = U[57] - U[63]; buffer27 = U[59] - U[61];
  204.     buffer18 = U[41] - U[47]; buffer28 = U[43] - U[45];
  205.     temp1 = buffer11 - buffer21;
  206.     temp2 = buffer12 - buffer22;
  207.     temp3 = buffer13 - buffer23;
  208.     temp4 = buffer14 - buffer24;
  209.     temp5 = buffer15 - buffer25;
  210.     temp6 = buffer16 - buffer26;
  211.     temp7 = buffer17 - buffer27;
  212.     temp8 = buffer18 - buffer28;
  213.     out[32] = temp1;
  214.     out[33] = temp2;
  215.     out[34] = temp3 - temp4;
  216.     out[35] = temp3 + temp4;
  217.     dummy1 = temp5 - temp7;
  218.     dummy2 = temp6 - temp8;
  219.     out[36] = dummy1 - dummy2;
  220.     out[37] = dummy1 + dummy2;
  221.     out[38] = -temp5 - temp7;
  222.     out[39] = temp6 + temp8;
  223.     temp1 = buffer11 + buffer21;
  224.     temp2 = buffer12 + buffer22;
  225.     temp3 = buffer13 + buffer23;
  226.     temp4 = buffer14 + buffer24;
  227.     temp5 = buffer15 + buffer25;
  228.     temp6 = buffer16 + buffer26;
  229.     temp7 = buffer17 + buffer27;
  230.     temp8 = buffer18 + buffer28;
  231.     out[40] = temp1;
  232.     out[41] = temp2;
  233.     out[42] = temp3 - temp4;
  234.     out[43] = temp3 + temp4;
  235.     dummy1 = temp5 - temp7;
  236.     dummy2 = temp6 - temp8;
  237.     out[44] = dummy1 - dummy2;
  238.     out[45] = dummy1 + dummy2;
  239.     out[46] = -temp5 - temp7;
  240.     out[47] = temp6 + temp8; /* 48 additions */
  241.  
  242.     temp1 = -U[1] - U[7];
  243.     temp2 = -U[33] - U[39];
  244.     temp3 = -U[17] - U[23];
  245.     temp4 = -U[49] - U[55];
  246.     temp5 = -U[9] - U[15];
  247.     temp6 = -U[25] - U[31];
  248.     temp7 = -U[57] - U[63];
  249.     temp8 = -U[41] - U[47];
  250.     out[48] = temp1;
  251.     out[49] = temp2;
  252.     out[50] = temp3 - temp4;
  253.     out[51] = temp3 + temp4;
  254.     dummy1 = temp5 - temp7;
  255.     dummy2 = temp6 - temp8;
  256.     out[52] = dummy1 - dummy2;
  257.     out[53] = dummy1 + dummy2;
  258.     out[54] = -temp5 - temp7;
  259.     out[55] = temp6 + temp8; /* 16 additions */
  260.  
  261.     temp1 = U[3] + U[5];
  262.     temp2 = U[35] + U[37];
  263.     temp3 = U[19] + U[21];
  264.     temp4 = U[51] + U[53];
  265.     temp5 = U[11] + U[13];
  266.     temp6 = U[27] + U[29];
  267.     temp7 = U[59] + U[61];
  268.     temp8 = U[43] + U[45];
  269.     out[56] = temp1;
  270.     out[57] = temp2;
  271.     out[58] = temp3 - temp4;
  272.     out[59] = temp3 + temp4;
  273.     dummy1 = temp5 - temp7;
  274.     dummy2 = temp6 - temp8;
  275.     out[60] = dummy1 - dummy2;
  276.     out[61] = dummy1 + dummy2;
  277.     out[62] = -temp5 - temp7;
  278.     out[63] = temp6 + temp8; /* 16 additions */
  279.     /*for (i=0; i<64; i++) printf("out[%d] = %f\n",i,out[i]);*/
  280.  
  281.     /********* Core Processing *********/
  282.     out1[0] = out[0];
  283.     out1[1] = out[1];
  284.     out1[2] = out[2];
  285.     out1[3] = 0.707106781 * out[3];
  286.     out1[4] = out[4];
  287.     out1[5] = 0.707106781 * out[5];
  288.     dummy1 = 1.306562964 * out[6];
  289.     dummy2 = 0.923879532 * (out[6] + out[7]);
  290.     dummy3 = -0.5411961 * out[7];
  291.     out1[6] = dummy1 - dummy2;
  292.     out1[7] = dummy2 + dummy3;   /* 5 multiplications, 3 additions */
  293.  
  294.     out1[8] = out[8];
  295.     out1[9] = out[9];
  296.     out1[10] = out[10];
  297.     out1[11] = 0.707106781 * out[11];
  298.     out1[12] = out[12];
  299.     out1[13] = 0.707106781 * out[13];
  300.     dummy1 = 1.306562964 * out[14];
  301.     dummy2 = 0.923879532 * (out[14] + out[15]);
  302.     dummy3 = -0.5411961 * out[15];
  303.     out1[14] = dummy1 - dummy2;
  304.     out1[15] = dummy2 + dummy3;   /* 5 multiplications, 3 additions */
  305.  
  306.     out1[16] = out[16];
  307.     out1[17] = out[17];
  308.     out1[18] = out[18];
  309.     out1[19] = 0.707106781 * out[19];
  310.     out1[20] = out[20];
  311.     out1[21] = 0.707106781 * out[21];
  312.     dummy1 = 1.306562964 * out[22];
  313.     dummy2 = 0.923879532 * (out[22] + out[23]);
  314.     dummy3 = -0.5411961 * out[23];
  315.     out1[22] = dummy1 - dummy2;
  316.     out1[23] = dummy2 + dummy3;   /* 5 multiplications, 3 additions */
  317.  
  318.     out1[32] = out[32];
  319.     out1[33] = out[33];
  320.     out1[34] = out[34];
  321.     out1[35] = 0.707106781 * out[35];
  322.     out1[36] = out[36];
  323.     out1[37] = 0.707106781 * out[37];
  324.     dummy1 = 1.306562964 * out[38];
  325.     dummy2 = 0.923879532 * (out[38] + out[39]);
  326.     dummy3 = -0.5411961 * out[39];
  327.     out1[38] = dummy1 - dummy2;
  328.     out1[39] = dummy2 + dummy3;   /* 5 multiplications, 3 additions */
  329.  
  330.     out1[24] = 0.707106781 * out[24];
  331.     out1[25] = 0.707106781 * out[25];
  332.     out1[26] = 0.707106781 * out[26];
  333.     out1[27] = 0.5 * out[27];
  334.     out1[28] = 0.707106781 * out[28];
  335.     out1[29] = 0.5 * out[29];
  336.     dummy1 = 0.923879532 * out[30];
  337.     dummy2 = 0.653281481 * (out[30] + out[31]);
  338.     dummy3 = -0.382683432 * out[31];
  339.     out1[30] = dummy1 - dummy2;
  340.     out1[31] = dummy2 + dummy3;   /* 9 multiplications, 3 additions */
  341.  
  342.     out1[40] = 0.707106781 * out[40];
  343.     out1[41] = 0.707106781 * out[41];
  344.     out1[42] = 0.707106781 * out[42];
  345.     out1[43] = 0.5 * out[43];
  346.     out1[44] = 0.707106781 * out[44];
  347.     out1[45] = 0.5 * out[45];
  348.     dummy1 = 0.923879532 * out[46];
  349.     dummy2 = 0.653281481 * (out[46] + out[47]);
  350.     dummy3 = -0.382683432 * out[47];
  351.     out1[46] = dummy1 - dummy2;
  352.     out1[47] = dummy2 + dummy3;   /* 9 multiplications, 3 additions */
  353.  
  354.     dummy1 = 1.306562964 * out[48];
  355.     dummy2 = 0.923879532 * (out[48] + out[56]);
  356.     dummy3 = -0.5411961 * out[56];
  357.     out1[48] = dummy1 - dummy2;
  358.     out1[56] = dummy2 + dummy3;
  359.     dummy1 = 1.306562964 * out[49];
  360.     dummy2 = 0.923879532 * (out[49] + out[57]);
  361.     dummy3 = -0.5411961 * out[57];
  362.     out1[49] = dummy1 - dummy2;
  363.     out1[57] = dummy2 + dummy3;  
  364.     dummy1 = 1.306562964 * out[50];
  365.     dummy2 = 0.923879532 * (out[50] + out[58]);
  366.     dummy3 = -0.5411961 * out[58];
  367.     out1[50] = dummy1 - dummy2;
  368.     out1[58] = dummy2 + dummy3;  
  369.     dummy1 = 1.306562964 * out[52];
  370.     dummy2 = 0.923879532 * (out[52] + out[60]);
  371.     dummy3 = -0.5411961 * out[60];
  372.     out1[52] = dummy1 - dummy2;
  373.     out1[60] = dummy2 + dummy3;  
  374.     dummy1 = 0.923879532 * out[51];
  375.     dummy2 = 0.653281481 * (out[51] + out[59]);
  376.     dummy3 = -0.382683432 * out[59];
  377.     out1[51] = dummy1 - dummy2;
  378.     out1[59] = dummy2 + dummy3;  
  379.     dummy1 = 0.923879532 * out[53];
  380.     dummy2 = 0.653281481 * (out[53] + out[61]);
  381.     dummy3 = -0.382683432 * out[61];
  382.     out1[53] = dummy1 - dummy2;
  383.     out1[61] = dummy2 + dummy3;  
  384.     temp1 = 0.5 * (out[54] + out[63]);
  385.     temp2 = 0.5 * (out[55] - out[62]);
  386.     temp3 = out[54] - out[63];
  387.     temp4 = out[55] + out[62];
  388.     temp5 = 0.35355339 * (temp3 - temp4);
  389.     temp6 = 0.35355339 * (temp3 + temp4);
  390.     out1[54] = temp1 - temp6;
  391.     out1[55] = temp2 + temp5;
  392.     out1[62] = temp5 - temp2;
  393.     out1[63] = temp1 + temp6; /* 22 multiplications 28 additions */
  394.     /*for (i=0; i<64; i++) printf("out1[%d] = %f\n",i,out1[i]);*/
  395.  
  396.     /********* Post Additions *********/
  397.     temp1 = out1[0] + out1[8];
  398.     temp2 = out1[1] + out1[9];
  399.     temp3 = out1[2] + out1[10];
  400.     temp4 = out1[3] + out1[11];
  401.     temp5 = out1[4] + out1[12];
  402.     temp6 = out1[5] + out1[13];
  403.     temp7 = out1[6] + out1[14];
  404.     temp8 = out1[7] + out1[15];
  405.     out[0] = temp1 + temp2;
  406.     out[1] = temp1 - temp2;
  407.     out[2] = temp4;
  408.     out[3] = temp3 - temp4;
  409.     out[4] = temp7 - temp6;
  410.     out[5] = temp8;
  411.     out[6] = -temp5 - temp7;
  412.     out[7] = temp6 + temp8;
  413.     temp1 = out1[0] - out1[8];
  414.     temp2 = out1[1] - out1[9];
  415.     temp3 = out1[2] - out1[10];
  416.     temp4 = out1[3] - out1[11];
  417.     temp5 = out1[4] - out1[12];
  418.     temp6 = out1[5] - out1[13];
  419.     temp7 = out1[6] - out1[14];
  420.     temp8 = out1[7] - out1[15];
  421.     out[8] = temp1 + temp2;
  422.     out[9] = temp1 - temp2;
  423.     out[10] = temp4;
  424.     out[11] = temp3 - temp4;
  425.     out[12] = temp7 - temp6;
  426.     out[13] = temp8;
  427.     out[14] = -temp5 - temp7;
  428.     out[15] = temp6 + temp8;
  429.     out[16] = out1[24] + out1[25];
  430.     out[17] = out1[24] - out1[25];
  431.     out[18] = out1[27];
  432.     out[19] = out1[26] - out1[27];
  433.     out[20] = out1[30] - out1[29];
  434.     out[21] = out1[31];
  435.     out[22] = -out1[28] - out1[30];
  436.     out[23] = out1[29] + out1[31];
  437.     temp1 = out1[16] - out1[24];
  438.     temp2 = out1[17] - out1[25];
  439.     temp3 = out1[18] - out1[26];
  440.     temp4 = out1[19] - out1[27];
  441.     temp5 = out1[20] - out1[28];
  442.     temp6 = out1[21] - out1[29];
  443.     temp7 = out1[22] - out1[30];
  444.     temp8 = out1[23] - out1[31];
  445.     out[24] = temp1 + temp2;
  446.     out[25] = temp1 - temp2;
  447.     out[26] = temp4;
  448.     out[27] = temp3 - temp4;
  449.     out[28] = temp7 - temp6;
  450.     out[29] = temp8;
  451.     out[30] = -temp5 - temp7;
  452.     out[31] = temp6 + temp8;
  453.     temp1 = out1[48] - out1[40];
  454.     temp2 = out1[49] - out1[41];
  455.     temp3 = out1[50] - out1[42];
  456.     temp4 = out1[51] - out1[43];
  457.     temp5 = out1[52] - out1[44];
  458.     temp6 = out1[53] - out1[45];
  459.     temp7 = out1[54] - out1[46];
  460.     temp8 = out1[55] - out1[47];
  461.     out[32] = temp1 + temp2;
  462.     out[33] = temp1 - temp2;
  463.     out[34] = temp4;
  464.     out[35] = temp3 - temp4;
  465.     out[36] = temp7 - temp6;
  466.     out[37] = temp8;
  467.     out[38] = -temp5 - temp7;
  468.     out[39] = temp6 + temp8;
  469.     out[40] = out1[56] + out1[57];
  470.     out[41] = out1[56] - out1[57];
  471.     out[42] = out1[59];
  472.     out[43] = out1[58] - out1[59];
  473.     out[44] = out1[62] - out1[61];
  474.     out[45] = out1[63];
  475.     out[46] = -out1[60] - out1[62];
  476.     out[47] = out1[63] + out1[61];
  477.     temp1 = -out1[32] - out1[48];
  478.     temp2 = -out1[33] - out1[49];
  479.     temp3 = -out1[34] - out1[50];
  480.     temp4 = -out1[35] - out1[51];
  481.     temp5 = -out1[36] - out1[52];
  482.     temp6 = -out1[37] - out1[53];
  483.     temp7 = -out1[38] - out1[54];
  484.     temp8 = -out1[39] - out1[55];
  485.     out[48] = temp1 + temp2;
  486.     out[49] = temp1 - temp2;
  487.     out[50] = temp4;
  488.     out[51] = temp3 - temp4;
  489.     out[52] = temp7 - temp6;
  490.     out[53] = temp8;
  491.     out[54] = -temp5 - temp7;
  492.     out[55] = temp6 + temp8;
  493.     temp1 = out1[40] + out1[56];
  494.     temp2 = out1[41] + out1[57];
  495.     temp3 = out1[42] + out1[58];
  496.     temp4 = out1[43] + out1[59];
  497.     temp5 = out1[44] + out1[60];
  498.     temp6 = out1[45] + out1[61];
  499.     temp7 = out1[46] + out1[62];
  500.     temp8 = out1[47] + out1[63];
  501.     out[56] = temp1 + temp2;
  502.     out[57] = temp1 - temp2;
  503.     out[58] = temp4;
  504.     out[59] = temp3 - temp4;
  505.     out[60] = temp7 - temp6;
  506.     out[61] = temp8;
  507.     out[62] = -temp5 - temp7;
  508.     out[63] = temp6 + temp8; /* 96 additions */
  509.     /*for (i=0; i<64; i++) printf("out[%d] = %f\n",i,out[i]);*/
  510.  
  511.     temp1 = out[0] + out[16];
  512.     temp2 = out[1] + out[17];
  513.     temp3 = out[2] + out[18];
  514.     temp4 = out[3] + out[19];
  515.     temp5 = out[4] + out[20];
  516.     temp6 = out[5] + out[21];
  517.     temp7 = out[6] + out[22];
  518.     temp8 = out[7] + out[23];
  519.     out1[0] = temp1 + temp3;
  520.     out1[1] = temp2 + temp4;
  521.     out1[2] = temp2 - temp4;
  522.     out1[3] = temp1 - temp3;
  523.     out1[4] = temp5;
  524.     out1[5] = temp6;
  525.     out1[6] = temp7;
  526.     out1[7] = temp8;
  527.     temp1 = out[8] + out[24];
  528.     temp2 = out[9] + out[25];
  529.     temp3 = out[10] + out[26];
  530.     temp4 = out[11] + out[27];
  531.     temp5 = out[12] + out[28];
  532.     temp6 = out[13] + out[29];
  533.     temp7 = out[14] + out[30];
  534.     temp8 = out[15] + out[31];
  535.     out1[8] = temp1 + temp3;
  536.     out1[9] = temp2 + temp4;
  537.     out1[10] = temp2 - temp4;
  538.     out1[11] = temp1 - temp3;
  539.     out1[12] = temp5;
  540.     out1[13] = temp6;
  541.     out1[14] = temp7;
  542.     out1[15] = temp8;
  543.     temp1 = out[8] - out[24];
  544.     temp2 = out[9] - out[25];
  545.     temp3 = out[10] - out[26];
  546.     temp4 = out[11] - out[27];
  547.     temp5 = out[12] - out[28];
  548.     temp6 = out[13] - out[29];
  549.     temp7 = out[14] - out[30];
  550.     temp8 = out[15] - out[31];
  551.     out1[16] = temp1 + temp3;
  552.     out1[17] = temp2 + temp4;
  553.     out1[18] = temp2 - temp4;
  554.     out1[19] = temp1 - temp3;
  555.     out1[20] = temp5;
  556.     out1[21] = temp6;
  557.     out1[22] = temp7;
  558.     out1[23] = temp8;
  559.     temp1 = out[0] - out[16];
  560.     temp2 = out[1] - out[17];
  561.     temp3 = out[2] - out[18];
  562.     temp4 = out[3] - out[19];
  563.     temp5 = out[4] - out[20];
  564.     temp6 = out[5] - out[21];
  565.     temp7 = out[6] - out[22];
  566.     temp8 = out[7] - out[23];
  567.     out1[24] = temp1 + temp3;
  568.     out1[25] = temp2 + temp4;
  569.     out1[26] = temp2 - temp4;
  570.     out1[27] = temp1 - temp3;
  571.     out1[28] = temp5;
  572.     out1[29] = temp6;
  573.     out1[30] = temp7;
  574.     out1[31] = temp8;
  575.     out1[32] = out[32] + out[34];
  576.     out1[33] = out[33] + out[35];
  577.     out1[34] = out[33] - out[35];
  578.     out1[35] = out[32] - out[34];
  579.     out1[36] = out[36];
  580.     out1[37] = out[37];
  581.     out1[38] = out[38];
  582.     out1[39] = out[39];
  583.     out1[40] = out[40] + out[42];
  584.     out1[41] = out[41] + out[43];
  585.     out1[42] = out[41] - out[43];
  586.     out1[43] = out[40] - out[42];
  587.     out1[44] = out[44];
  588.     out1[45] = out[45];
  589.     out1[46] = out[46];
  590.     out1[47] = out[47];
  591.     out1[48] = out[48] + out[50];
  592.     out1[49] = out[49] + out[51];
  593.     out1[50] = out[49] - out[51];
  594.     out1[51] = out[48] - out[50];
  595.     out1[52] = out[52];
  596.     out1[53] = out[53];
  597.     out1[54] = out[54];
  598.     out1[55] = out[55];
  599.     out1[56] = out[56] + out[58];
  600.     out1[57] = out[57] + out[59];
  601.     out1[58] = out[57] - out[59];
  602.     out1[59] = out[56] - out[58];
  603.     out1[60] = out[60];
  604.     out1[61] = out[61];
  605.     out1[62] = out[62];
  606.     out1[63] = out[63]; /* 64 additions */
  607.     /*for (i=0; i<64; i++) printf("out1[%d] = %f\n",i,out1[i]);*/
  608.  
  609.     temp1 = out1[0] + out1[32];
  610.     temp2 = out1[1] + out1[33];
  611.     temp3 = out1[2] + out1[34];
  612.     temp4 = out1[3] + out1[35];
  613.     temp5 = out1[4] + out1[36];
  614.     temp6 = out1[5] + out1[37];
  615.     temp7 = out1[6] + out1[38];
  616.     temp8 = out1[7] + out1[39];
  617.     U[0] = temp1 + temp5;
  618.     U[8] = temp2 + temp6;
  619.     U[16] = temp3 + temp7;
  620.     U[24] = temp4 + temp8;
  621.     U[32] = temp4 - temp8;
  622.     U[40] = temp3 - temp7;
  623.     U[48] = temp2 - temp6;
  624.     U[56] = temp1 - temp5;
  625.     temp1 = out1[8] + out1[40];
  626.     temp2 = out1[9] + out1[41];
  627.     temp3 = out1[10] + out1[42];
  628.     temp4 = out1[11] + out1[43];
  629.     temp5 = out1[12] + out1[44];
  630.     temp6 = out1[13] + out1[45];
  631.     temp7 = out1[14] + out1[46];
  632.     temp8 = out1[15] + out1[47];
  633.     U[1] = temp1 + temp5;
  634.     U[9] = temp2 + temp6;
  635.     U[17] = temp3 + temp7;
  636.     U[25] = temp4 + temp8;
  637.     U[33] = temp4 - temp8;
  638.     U[41] = temp3 - temp7;
  639.     U[49] = temp2 - temp6;
  640.     U[57] = temp1 - temp5;
  641.     temp1 = out1[16] + out1[48];
  642.     temp2 = out1[17] + out1[49];
  643.     temp3 = out1[18] + out1[50];
  644.     temp4 = out1[19] + out1[51];
  645.     temp5 = out1[20] + out1[52];
  646.     temp6 = out1[21] + out1[53];
  647.     temp7 = out1[22] + out1[54];
  648.     temp8 = out1[23] + out1[55];
  649.     U[2] = temp1 + temp5;
  650.     U[10] = temp2 + temp6;
  651.     U[18] = temp3 + temp7;
  652.     U[26] = temp4 + temp8;
  653.     U[34] = temp4 - temp8;
  654.     U[42] = temp3 - temp7;
  655.     U[50] = temp2 - temp6;
  656.     U[58] = temp1 - temp5;
  657.     temp1 = out1[24] + out1[56];
  658.     temp2 = out1[25] + out1[57];
  659.     temp3 = out1[26] + out1[58];
  660.     temp4 = out1[27] + out1[59];
  661.     temp5 = out1[28] + out1[60];
  662.     temp6 = out1[29] + out1[61];
  663.     temp7 = out1[30] + out1[62];
  664.     temp8 = out1[31] + out1[63];
  665.     U[3] = temp1 + temp5;
  666.     U[11] = temp2 + temp6;
  667.     U[19] = temp3 + temp7;
  668.     U[27] = temp4 + temp8;
  669.     U[35] = temp4 - temp8;
  670.     U[43] = temp3 - temp7;
  671.     U[51] = temp2 - temp6;
  672.     U[59] = temp1 - temp5;
  673.     temp1 = out1[24] - out1[56];
  674.     temp2 = out1[25] - out1[57];
  675.     temp3 = out1[26] - out1[58];
  676.     temp4 = out1[27] - out1[59];
  677.     temp5 = out1[28] - out1[60];
  678.     temp6 = out1[29] - out1[61];
  679.     temp7 = out1[30] - out1[62];
  680.     temp8 = out1[31] - out1[63];
  681.     U[4] = temp1 + temp5;
  682.     U[12] = temp2 + temp6;
  683.     U[20] = temp3 + temp7;
  684.     U[28] = temp4 + temp8;
  685.     U[36] = temp4 - temp8;
  686.     U[44] = temp3 - temp7;
  687.     U[52] = temp2 - temp6;
  688.     U[60] = temp1 - temp5;
  689.     temp1 = out1[16] - out1[48];
  690.     temp2 = out1[17] - out1[49];
  691.     temp3 = out1[18] - out1[50];
  692.     temp4 = out1[19] - out1[51];
  693.     temp5 = out1[20] - out1[52];
  694.     temp6 = out1[21] - out1[53];
  695.     temp7 = out1[22] - out1[54];
  696.     temp8 = out1[23] - out1[55];
  697.     U[5] = temp1 + temp5;
  698.     U[13] = temp2 + temp6;
  699.     U[21] = temp3 + temp7;
  700.     U[29] = temp4 + temp8;
  701.     U[37] = temp4 - temp8;
  702.     U[45] = temp3 - temp7;
  703.     U[53] = temp2 - temp6;
  704.     U[61] = temp1 - temp5;
  705.     temp1 = out1[8] - out1[40];
  706.     temp2 = out1[9] - out1[41];
  707.     temp3 = out1[10] - out1[42];
  708.     temp4 = out1[11] - out1[43];
  709.     temp5 = out1[12] - out1[44];
  710.     temp6 = out1[13] - out1[45];
  711.     temp7 = out1[14] - out1[46];
  712.     temp8 = out1[15] - out1[47];
  713.     U[6] = temp1 + temp5;
  714.     U[14] = temp2 + temp6;
  715.     U[22] = temp3 + temp7;
  716.     U[30] = temp4 + temp8;
  717.     U[38] = temp4 - temp8;
  718.     U[46] = temp3 - temp7;
  719.     U[54] = temp2 - temp6;
  720.     U[62] = temp1 - temp5;
  721.     temp1 = out1[0] - out1[32];
  722.     temp2 = out1[1] - out1[33];
  723.     temp3 = out1[2] - out1[34];
  724.     temp4 = out1[3] - out1[35];
  725.     temp5 = out1[4] - out1[36];
  726.     temp6 = out1[5] - out1[37];
  727.     temp7 = out1[6] - out1[38];
  728.     temp8 = out1[7] - out1[39];
  729.     U[7] = temp1 + temp5;
  730.     U[15] = temp2 + temp6;
  731.     U[23] = temp3 + temp7;
  732.     U[31] = temp4 + temp8;
  733.     U[39] = temp4 - temp8;
  734.     U[47] = temp3 - temp7;
  735.     U[55] = temp2 - temp6;
  736.     U[63] = temp1 - temp5; /* 128 additions */
  737. }
  738.